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Abstract: A geometric derivation of numerical integrators for nonholonomic 
systems and optimal control problems is obtained. It is based in the classical 
technique of generating functions adapted to the special features of nonholonomic 
systems and optimal control problems. 
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1. INTRODUCTION 

Standard methods for simulating the motion of 
a dynamical system usually ignore many of the 
geometric features of this system (simplecticity, 
conservation laws, symmetries...). However, new 
methods have been recently developed, called geo- 
metric integrators, which are concerned with some 
of the extra features of geometric nature of the 
dynamical system (see [HaLuWa:02]). 

In the first part of the paper, we propose a class 
of geometric integrators for nonholonomic systems 
[Leomar:96D,NeiFuf:72] based on a discretization 
of the Lagrangian function (in a more precise 
sense, we discretize the action function) and a co- 
herent discretization of the constraint forces (see 
[LeMaSa:02]). These equations will be conceptu- 
ally equivalent to the proposed for systems with 
external forces (see [MarWes:01]). Finally, second 
part corcerns with the construction of symplectic 
integrators for optimal control theory by using 
generating functions of the second kind. 
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2. NONHOLONOMIC SYSTEMS 

2.1 Geometrical formulation of nonholonomic systems 

Let Q be a n-dimensional differentiable manifold, 
with local coordinates [q^) and tangent bundle 
TQ, with induced coordinates {q^q^). Consider a 
Lagrangian system, with Lagrangian L : TQ — > R, 
subject to nonholonomic constraints, defined by a 
submanifold D of the velocity phase space TQ. 
We will assume that dimZ? = 2n — m and that 
D is locally described by the vanishing of m inde- 
pendent functions 0° (the "constraint functions"), 

satisfying the rank condition rank — 

V 

In the sequel, we will follow a Hamiltonian point 
of view. The canonical coordinates on T*Q (the 
cotangent bundle of Q) are denoted by {q^,Pi). 
Assume, for simplicity, that the Lagrangian L 
is hyperregular, that is, the Legendre transfor- 
mation Leg : TQ T*Q,{q\q') ^ {q\p, ^ 
dL/dcf), is a global diffeomorphism. The con- 
straint functions on T*Q become ^f" — 4'"-oLeg~^ , 

i.e. '^°-{q\pi) = 0°((?', -^) , where the Hamilto- 
nian H : T*Q ^ R \s defined hy H ^ El o 
Leg~^. Here, El denotes the energy of the sys- 
tem, locally defined hy El = 'i^§§r — L. Since 



locally Leg ^{q^,Pi) 



dH 

{<l\^), then H =p,cf 



Liq^,q^) , where is expressed in terms of and 
Pi by using Leg~^. 

The equations of motion for the nonholonomic 

system on T*Q can now be written as follows (see 
[CaLeMa:99,Marl:95] and references therein) 



■ i dH 
1 = a- 

OPi 

dH ^ 9*" 

Pi = ~in ~ '^"^ — 

aq^ opj 



(1) 



together with the constraint equations \E'°(g,p) = 
0, where Hij are the components of the inverse of 
the matrix (W^) = {d^ H / dpidpj) . Note that 

— ^3i){q,P) = (-^ ° Leg )[q,p). 



dp 



dqi 



Let M denote the image of the constraint subman- 
ifold D under the Legendre transformation, and 
let F be the distribution on T*Q along M, whose 
annihilator is given by F° = Legt,{F°)). Here, F° 
represents the constraint forces subbundle, locally 
defined by 



F° = span{^" 



dcj)" 



dq'} 



The Hamiltonian equations of motion of the non- 
holonomic system can be then rewritten in intrin- 
sic form as 



{ixuJQ - dH)\M e F" 



(2) 



where ujq = —cLOq = dq^Adpi (with 9q = pi dq^) is 
the canonical symplectic form on T*Q. Suppose in 
addition that the following compatibility condition 
F-L n TM = {0} holds, where " _L " denotes the 
symplectic orthogonal with respect to u)q. Ob- 
serve that, locally, this condition means that the 

matrix (C"h = i^U,, ^— is regular. The 

\dPi / dpj J 
compatibility condition is not too restrictive, since 
it is trivially verified by the usual systems of me- 
chanical type (i.e. with a Lagrangian of the form 
kinetic minus potential energy), where the Tiij 
represent the components of a positive definite 
Riemannian metric. The compatibility condition 
guarantees, in particular, the existence of a unique 
solution of the constrained equations of motion (2) 
which, henceforth, will be denoted by Xh,m on 
the Hamiltonian side and Leg~^{XH,M) = ^l,d 
on the Lagrangian side. 

Moreover, if we denote by Xh the Hamiltonian 
function of H, i.e., Ixh'^Q = dH then, using 
the constraint functions, we may explicitely de- 
termine the Lagrange multipliers Aa as = 
-Ca6X/f(*'') . Next, writing the 1-form A = 
— CafcXjf (5''')^-7iji(i(7* then, the nonholonomic 
equations are equivalently rewritten as 



q = 
Pi 



dH 

dpi 
dH . 

— TT- ~ ■ 

dq^ 



(3) 



for initial conditions {qo,Po) G M and A = Aj dg*. 
We also denote by A = Leg* (A) the 1-form on 
TQ wich represents the constraint force once the 
Lagrange multipliers have been determined. 

Now, consider the flow Ft : M ^ M, t € 
/ C M of the vector field Xh.Mi solution of the 
nonholonomic problem. Since (3) is geometrically 
rewritten as 



ixH.M'^Q = dH + A, 



then 



Lxh.m^Q = diixn.M^Q -H)-A, 
or, equivalently, Lxh,m(^Q = d{L o Leg~'^) 



Therefore, integrating 
F'hSQ-OQ=d \ I LoFtdt 



Ft* A, 



A . 



(4) 



where Ft is the flow of the vector field ^l,d- 



2.2 "Generating functions" and nonholonomic 
mechanics 

In what follows, we will follow similar arguments 
for the construction of generating functions for 
symplectic or canonical maps [Arn;78]. However, 
because of equation (4), we have that the non- 
holonomic flow is not a canonical transformation; 



I.e., 



F^UIQ - WQ 




(5) 



This description will allow us to construct a new 
family of nonholonomic integrators for equations 
(19). Denote by tt^ : T*Q x T*Q T*Q, i = l,2, 
the canonical projections. Consider the following 
forms 



6: 



-dQ . 



Denote by ip^ ■ Graph(F,,) ^ T*Q x T*Q the 
inclusion map and observe that Graph (F/j) C Mx 
M. Then, from (4) i*p 6 is equal to 



(tti 



iGraph(Ffc) 



L o Ftdt 



Ft*A 



Let ((JojPOj giiPi) be coordinates in T*Q x T*Q 
in a neighborhood of some point in Graph(Fft,). If 
(90, Po, 91, Pi) e Graph(F;i) then *"(9o,Po) = 
and '^"■{qi,pi) = 0. Moreover, along Graph(F/j), 

qi = qi{qo,po), pi = pi{qo,po) and 



pidqi -pQdqo = d^J L{q{t),q{t)) 

- f l{q{t),q{t)), 
Jo 



dt 



(6) 



where {q{t),q{ty) = Ft{qo,qo) with Leg{qo,qo) = 
(90, Po)- Here, Ft denotes the flow of ^l,d- Equa- 
tion (6) is satisfied along Graph(F^). 

Assume that, in a neighborhood of some point 
X e Graph(F/(), we can change this system of 
coordinates to a new coordinates {qo,qi)- Denote 
by 

S''{qa,qi) = / L{q{t),q{t))dt, 
JQ 

where q{t) is a solution curve of the nonholonomic 
problem with q{0) = q and q{h) = q\ and an 
adequate extension of . It is easy to show that 
this solution always exists for adequate values of 
qo and qi. 



Thus, we deduce that 

Qgh r-h 



Po 



Pi = 



dqa 
dS^ _ 
dqi 



A{q{t),m)^ 
ctlo 

'A{q{t),q{t))^, 
oqi 



(7) 



where {qo,qi) verifies the constraint functions 
V^ilo, Qi,h) = 0, explicitely defined by 

ip°-{qo,qi,h) = 



tively, q{h) = qk and q{2h) = qk+i) for the first 
integral (resp., second integral) of the right-hand 
side. 

Proof: It is suffices to prove the result for A'' = 2; 
that is, 

S^^{qo,q2) = S^{q^,q,)+S''{quq2), 
where qi verifies condition (9). 
Since 

Pidqi-podqo = dS^{qo,qi) - / A{q{t),q{t)) , 

Jo 

Pi dqi - pi dqi = dS''{quqi) - / A(?(t), q{t)) , 

Jh 

then 

Pi dq2 -podqo = d{S''iqo,qi) + S^ {qi , 92 )) 

- / A{q{t),m)- / mt),m)- 

Jo Jh 

Since the variables qi do not appear on the left- 
hand side term, we obtain expression (9). More- 
over, for this choice of qi then S^^{qo,q2) = 
S'^iqoTqi) + S'^iqi^qi) is a "generating function 
of the first kind" of i=2/,. I 

Equations (9) determine an implicit system of 
difference equations which permit us to obtain q2 
from the initial data go and qi. 



dS f ~ dq 

*''(go,--^(go,gi) + / A(g(f),g(t)) — ), (8) 2.3 Nonholonomic integrators 
OQo Jo ^90 



where q{t) is a solution of the nonholonomic 
problem with g(0) = qo and q{h) = qu- 

Next, we will show how the group composite law 
of the flow Fh, Fj^fi = Fh o . . . o Fh, is expressed in 

V ' 

N 

terms of the corresponding "generating functions" 
S^. Moreover, the following Theorem will result 
in a new construction of numerical integrators 
for nonholonomic mechanics when we change the 
"generating function" and the constraint forces by 
appropriate approximations. 

Theorem 2.1. The function S'^^, the "generating 
function" for -Fjvh, is given by 

N-l 

fe=o 

where qk, 1 < k < N — 1, are points verifying 

D2S'\qk-i,qk) + D,S'\qk,qk+l) = 
£A{q{t),q{t))-l^+rHm,m)-l^, (9) 



In the sequel and, for simplicity, assume that Q is 
a vector space. Since S^{qo, qi) = Jq L{q{t), q{t)) dt, 
where q{t) is a nonholonomic solution with q{0) = 
qo and q{h) = qi, we can obtain nonholonomic 
integrators by taking adequate approximations of 
the "generating function" S''* and the extra-term 

J,^A{q{t),q{t))- 

Consider, for instance, the approximation 



Sa{<lo,qi) =hL{{l-a)qo + aqi, 



Qi-Qo . 
h ' 



(10) 



dq 



1 Jh 



dqo 



and q{t) is a solution curve of the nonholonomic 
problem with g(0) = qk-i and q{h) = qk (respec- 



for some parameter a G [0, 1]. (In general, we will 
write 55(go,gi)«5''(go, 91).) 

A natural approximation of the constraint forces 
adapted to our choice of approximation for S'^ are 



\iqit),qit))p- 
oqo 

(1 - a)/iA((l - a)qo + aqi, ^ ^° ) , 
A(g(t), g(t))^ ~ ahA{{l - a)qo + aq,, ^^) 



Consequently, we obtain the following numerical 
method for nonholonomic systems 



D2S^{qk-i,qk) + D,S^{qk,qk+l) ^ 

ulim \ I Qk-qk-is 
ahK{(l ~ a)qk-i + aqk, ) 

+ (1 - a)/iA((l - a)qk + aqu+i, *±1-— ^) ^ 



with 1 < fc < A^ — 1 and initial condition satisfying 



dqo 



(90,91) 



+ (1 - a)hAiil - a)qo + aq,, ^^-^)) = 



Example 2.2. Nonholonomic particle. 

Consider the Lagrangian L : TR'^ R 

subject to the constraint (j) — z — yx = 0. Taking 
a = 1/2 in (10) we obtain a geometric integrator 
for the continuous nonholonomic problem. The 
first figure compares the method introduced here 
to the traditional Runge-Kutta method of fourth 
order, showing an improvement in several orders 
of magnitude. Observe that, in this scale, the 
value of the energy in each step of our algorithm 
is practically undistinguishable from the initial 
value of the energy. 




450 500 



The second figure is a comparison between our 
method and the one proposed in [CorMar:01]. 
A similar behaviour is observed. Nevertheless, a 
slightly better behaviour can also be appreciated, 
where the proposed algorithm shows on average a 
better preservation of the original energy. 



3. OPTIMAL CONTROL THEORY 

3. 1 Geometric formulation of optimal control problems 

A general optimal control problem consists of a 
set of differential equations 

q' ^T\q{t),u{t)),l<i<n, (11) 

where denote the states and u the control 
variables, and a cost function L(g, u). Given some 
boundary conditions (usually go — qlto) and qp = 
q{tf)) the aim is to find a C^-piecewise smooth 
curve c{t) — {q{t),u{t)), satisfying the control 
equations (11) and minimizing the functional 

J{c)^ r L{q{t),u{t))dt. (12) 



In a global description, one assumes a fiber bundle 
structure tt : U — > Q, where Q is the configu- 
ration manifold with local coordinates q^ and U 
is the bundle of controls, with local coordinates 
{q\ u"), l<«<n, l<a<TO. 

The ordinary differential equations (11) on Q 
depending on the parameters u can be seen as 
a vector field F along the projection map tt, that 
is, F is a smooth map F : U — > TB such that the 
diagram 




is commutative. This vector field is locally written 

asr = r(<7,u)A. 

The solutions of such problem are provided by 
Pontryaguin's maximum principle. If we construct 
the Hamiltonian function 

H{q,p,u)^L{q,u)+p,r{q,u) (13) 

where Pi, 1 < i < n, are now considered as 
Lagrange's multipliers, then a curve 7 ; R ^ [/, 
7(t) = {q{t),u{t)) is an optimal trajectory if there 
exists functions Pi{t), 1 < i < n such that they 
are solutions of the Hamilton equations: 



q\t) = —{q{t).P{t)Mt)) 
dH 

Mt)^-g^{q{t),p(.t),uit)) 



(14) 



and 



H{q{t),p{t),u{t))^ mm H{q{t),p{t),v), (15) 

V 

with t G [tQ,tf]. This last condition is usually 
replaced by 

'"^ - - - (16, 



— 0, 1 < a < m 



when we are looking for extremal trajectories. 



It is well known that the Pontryaguin's neces- 
sary conditions for cxtrcmality liavc a geometric 
interpretation in terms of presymplectic hamilto- 
nian system. The total space of the system will 
be T*Q Xq U. Let luq be the canonical sym- 
plectic form on T*Q and consider the canonical 
projection pr^ : T*Q Xq U — > T*Q. Denote 
by OJ = prlujQ the induced closed 2-form on 
T*Q Xq U. The 2-form lu is degenerate and its 
characteristic distribution is locally spanned by 
d/du'^, 1 < a < m. Define the Pontryaguin's 
hamiltonian function H : T*Q xq U — > M as 



follows H{aq 



= L(u„ 



-aq{r{uq)) where a, € 



T*Q and Uq e pr (q). Obviously, the coordinate 
expression of H is (13). 

Equations (14) (15) and (16) are intrinsically 
written as 



ix(jO = dH 



(17) 



Applying the Dirac-Bergmann-Gotay-Nester al- 
gorithm to the presymplectic system {T*Q Xq 
U, Q, H) we obtain that equations (16) correspond 
to the primary constraints for the presymplec- 
tic system: 6°- = = 0. The equations have 
solution along the first constraint submanifold 
Pq determined by the vanishing of the primary 
constraints. On the points of Pq there is at least 
a pointwise solution of Equation (17) but such 
solutions are not, in general, tangent to Pq. These 
points must be removed leaving a subset Pi C Pq 
(it is assumed tan Pi is also a submanifold). Then, 
we have to restrict Pi to a submanifold where the 
solutions of (17) are tangent to Pi. Proceeding 
further, we obtain a sequence of submanifolds 



Pk 



P2^Pi^Po^T*QxqU 



If this algorithm stabilizes, i.e. there exists a 
positive integer /c G N such that Pfe = P^+i 
and dimPk ^ 0, then we will obtain an stable 
submanifold Pf = Pk, on which a vector field 
exists such that 



(18) 



The constraints determining Pj are known in the 
control literature as higher order conditions for 
optimality. Therefore, a necessary condition for 
optimality of the curve; 7 : R ^ [/, ^{f?) = 
{q{t),u{t)) will be the existence of a lift 7 of 7 
to Pf such that 7 will be an integral curve of a 
solution of Equations (18). 

In the regular case, the final constraint algo- 
rithm is Pq (that is, Pq = P/) and all the con- 
straints are second class following the classical 
classification of Dirac. In such case (Po,wo) is a 
symplectic manifold, where f2o denotes the re- 
striction of the presymplectic 2-form to the con- 
straint submanifold Pq. Locally, the symplecticity 
of (PcCJo) is equivalent to the regularity of the 



matrix 



. The dynamical equa- 



du'^dvP / 1 ^„ . ^„ 

^ / l<.a,&<.m 

tions for the optimal control problem will be 

ixp^^Q = dH\p^ (19) 

Taking coordinates {q^,Pi) on Pq, then the dynam- 
ical equations are: 




dHi. 



dq 



(20) 



where we have substituted in (14) the control 
variables by its value u'^ = f°'{q,p) apply- 
ing the implicit function theorem to the primary 
constraints (/)° = 0. In such case, there exists a 
unique solution Xp^ of Equation (19) and its flow 
preserves the symplectic 2-form wq, i.e. it is a 
canonical transformation. 



3.2 Generating functions of the second kind 

Let {A4, n) be an exact symplectic manifold (fl is 
symplectic and exact, = — dO) and suppose that 
P : — > is a transformation from Ai to itself 
and Graph(P) the graph of P, Graph(P) C M x 
M. Denote hy ni : M x M M, i = 1,2 the 
canonical projections and the forms: 



6: 



:7r*e -7r*e 



= 7r217-7rjli: 



-de 



Denote hy ip ■ Graph(P) ^ AixAi the inclusion 
map. Then, P is a canonical transformation if 
and only if i*pCl = 0, that is, if Graph(P) is a 
lagrangian submanifold of {M x M,^). In such a 
case, i*pQ, = —dipQ = and, at least locally, there 
exists a function S : Graph P ^ K such that 

i*pe = dS (21) 

Taking {q^,Pi) as natural coordinates in Graph(P) 
and {q^,Pi,<t,Pi) the coordinates in M x A4, 
then, along Graph(P), q* = cC{q,p) and pi = 
Pi(9)P) and Pi dq* — Pidq^ = dS{q,p). Suppose 
that (g\Pi) are independent local coordinates on 
Graph(P) (see [Arn:78]); i.e. S = S{q,p) Since 

Pi dq" - Pi dq"- = -q' dp^ + d(q'pi) - Pi dq^ = dS, 

if we define S2{q,p) = q*Pi — S{q,p), where p 
is expressed in terms of p and q, then q* dpi + 
Pidq^ = dS2{q,p) 

Definition 3.1. The function S2{q, p) will be called 
a generating function of the second kind of 
the canonical transformation P. 

Now, suppose that {M.,Q,,H) is a hamiltonian 

system and Xh its hamiltonian vector field, say 
ixh^ = dH. Denote hy Fhi M ^ M its flow. 



Theorem 3.2. Let a function be defined by 



N-l 



k=0 

where qk, ^ < k < N, and pk, < k < N — 1, are 
stationary points of the right-hand side, that is 

qk+i = -g^{qk,Pk+i), o<fc<iV-i 

Pk = -g^{qk,Pk+i), 0<k<N-l 



then is a generating function of the second 
kind for F^h : M ^ M. 

Proof: It is similar to that of Theorem 2.1. | 
Finally, we have the following 

Proposition 3. 3. A generating function of the sec- 
ond kind for is given by 

S2{qo,Ph)=Phqh- I (pdq-Hdt) 
Jo 

where t — > {q{t),p{t)) is an integral curve of 
the Hamilton equations such that g'(O) = go and 
p{h) =ph. 



3.3 Generating functions of the second kind and 
discrete optimal control problems 

Prom Proposition 3.3 a generating function of 

the second kind for the Hamiltonian system 
(Po, fio, -ff|Po) which determines the dynamics of 
the optimal control problem given by (11) and 
(12) is 

S2{qQ,Ph) = PhQh 

- f\p{t)q{t)-H\pMi)^P{t))) dt (22) 
Jo 

where t {q{t),p{t)) is an integral curve of the 
vector field Xp„ with (g(0),p(0)) = {qo,Po) and 

{q{h),p{h)) = {qh,Ph)- 

We now turn to the construction of a numerical in- 
tegrator for the Hamiltonian system (Pqj ^Oi H\p^^) 
by using an approximation of the generating func- 
tion. The proposed methods also realize the inte- 
gration steps by canonical transformations; there- 
fore, they are symplectic integrators. 

Example 3.4- Consider, for instance, the following 
approximation to 82. 



'S'2 (i?fc,Pfc+i) = Pk+iqk+i - hp, 



k+l 



qk+1 — qk 



h 



where Ld and Td are adequate approximations to 
i|P(, and Fipq, respectively. 

Denote by f{qk,Pk+i) the function f{qk,Pk+i) = 

hTd{qk,Pk+i)+qk- Since ^''^^^ = ^d{qk,Pk+i) 
then, 

S2{qk,Pk+i) = Ld{qk,Pk+i) +Pk+if{qk,Pk+i) 
and hence the equations 

,qk,Pk+i) 

(23) 



Pk = 



qk+1 



dq 

dp 



{qk,Pk+i) 



+hLd{qk,Pk+i) + hpk+iTd{qk,Pk+i) 



are exactly the discrete equations corresponding 
to the classical discrete optimal control problem 
(see [Lew:86]), determined by the control equa- 
tions: qlj^^ = r{qk,Uk), ((go) given) and with 

associate perfomance index: J = J2k^o ^<i{qk, Uk) 
Observe that this discrete optimal control prob- 
lem is symplectic in the sense explained in the 
subsection above. 
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